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Abstract 

Singular spectrum analysis is developed as a nonparametric spec¬ 
tral decomposition of a time series. It can be easily extended to the de¬ 
composition of multidimensional lattice-like data through the filtering 
interpretation. In this viewpoint, the singular spectrum analysis can 
be understood as the adaptive and optimal generation of the filters and 
their two-step point-symmetric operation to the original data. In this 
paper, we point out that, when applied to the multidimensional data, 
the adaptively generated filters exhibit symmetry properties result¬ 
ing from the bisymmetric nature of the lag-covariance matrices. The 
eigenvectors of the lag-covariance matrix are either symmetric or an¬ 
tisymmetric, and for the 2D image data, these lead to the differential- 
type filters with even- or odd-order derivatives. The dominant filter 
is a smoothing filter, reflecting the dominance of low-frequency com¬ 
ponents of the photo images. The others are the edge-enhancement 
or the noise filters corresponding to the band-pass or the high-pass 
filters. The implication of the decomposition to the image denoising 
is briefly discussed. 
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1 Introduction 


Singular spectrum analysis (SSA) is a nonparametric and adaptive decom¬ 
position of a time series. In the SSA, the original time series is decomposed 
into arbitrary number of additive components and they are interpreted as the 
slowly varying trend, oscillatory and noise components, respectively [2,3]. In 
the conventional approach, the original time series { xf, i = 0,..., L — 1} is em¬ 
bedded into the trajectory matrix X through the lagged vectors of length K 
by forming L — K + 1 lagged vectors. Then the singular value decomposition 
of the trajectory matrix X is carried out as 


X = Ad + X 2 + .... + X K , 


( 1 ) 


with Xj are the rank 1 matrices. After grouping the matrices A j, the time se¬ 
ries is reconstructed by projecting the matrices Xj into Hankel form through 
the anti-diagonal averaging. Recently, it has been shown that these SSA 
procedure can equivalently be interpreted as the two-step (forward and the 
reverse) filtering of a time series [6,8]. The eigenvector {v^;k = 1 
of the matrix X T X can be interpreted as the convolution coefficients of the 
non-causal FIR filter. Let us denote the discrete Fourier transform of the 
eigenvector for X T X as {v^}. In Fourier space, the first (second) 

step filtering can be expressed as the multiplication of v^* ( v to the 
original sequences. Then the two-step filtering can be expressed as the mul¬ 
tiplication of the zero-phase Liters \v^\ 2 [7]. The normalization and the 
perfect-reconstruction properties can be expresses as 
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and these are the consequences of the normalization and the completeness 
properties of the eigenvectors {u^} for X T X [8]. 

Though the filtering interpretation of the SSA algorithm is equivalent to 
the conventional treatment, it has several advantages. In the conventional 
approach, the singular value decomposition of the trajectory matrix and the 
reconstruction procedure through the Hankel projection seems to be quite 
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different procedures and this makes it difficult to extend SSA to higher di¬ 
mension. With the filtering interpretation, these procedures can be seen 
as the forward and the reverse two-step point-symmetric convolution with 
respect to the reference points. This largely simplify the algorithms and en¬ 
ables us to extend the SSA to the decomposition of multidimensional data 
with arbitrary dimension [9]. Also, we can study the properties of the filters 
separately and it helps us for the detailed study of the spectral meaning of 
the decomposition. 

There have been several works in which the SSA algorithm has been 
applied to 2D image data [4,12] or three dimensional polygonal data [10,11]. 
In this paper, we focus on the SSA decomposition of the 2D image data. 
Under the periodic boundary condition for the image data, the matrix X T X 
becomes to be the bisymmetric type and it leads to either the symmetric or 
the antisymmetric eigenvectors. A single step filtering generated with these 
eigenvectors can be interpreted as the differential operation to the images 
and the symmetric (antisymmetric) eigenvalues lead to the derivative filters 
of the even- ( odd- ) order. In the image processing, various differential filters 
are designed to detect and enhance the edge lines in the images. In the SSA 
decomposition of the 2D image data, these differential filters are adaptively 
and optimally generated. We have examined these points in detail, focusing 
our attention to the characteristics of the filters in the SSA. Finally, we briefly 
discuss the possibility of the denoising with the SSA decomposition. 

In Sec. 2, we give a brief summary of the SSA algorithm in 2D image 
decomposition, emphasizing the role of filters. In Sec. 3, we give an example 
of the SSA decomposition of the image data, and the symmetry of the eigen¬ 
vector leads to the symmetry properties of the filters and these are reflected 
as the filter characteristics for the image data. An implication of these prop¬ 
erties to the image denoising is briefly discussed. In Sec. 4, we summarize 
the results of this paper. 

2 SSA Decomposition of 2D Image Data 

2.1 Basic algorithm 

Let us consider the 2D image data given by the matrix A = {a^} with 
0 < i < M — 1 and 0 < j < N— 1. We adopt the moving window of any shape. 
For simplicity we assume here the m x n rectangular form (m < M, n < N ), 
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which is shown as the box [• • •] in the following equation 
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from left to right and from top to bottom, for 
example, the trajectory matrix A" can be defined as 
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To introduce the filter for the 2D data, we consider an arbitrary K = mn 

dimensional orthonormal vector {vj : j = 1 = mn} . We arbitrarily 

assign the components of Vj in the moving window, 
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In the SSA, we adopt the vector set : k — 1,..., K = mn} which is the 
eigenvectors of the matrix X T X and the filters are defined as, 
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The following convolution-type linear relation express the forward filtering 

v A k ) 
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= v\ ’aij + v\ 'a i+ ij + ... + VK_ m+1 a itj+n -i 

T... T Vj} Oi+m—1 ,j+n —1 • (8) 

Next, we perform the reverse filtering with the point symmetric filter 
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for the reference point a^- with the relation 
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We assume the periodic boundary for the image data A. In the SSA, the 
vectors {v^} are determined as the eigenvector of the lag-covariance matrix 
C = X t X/(L — K + 1). It is easy to show that the completeness of the 
orthonormal eigenvectors leads to the perfect-reproduction property of the 
two-step filtering irrespective of the window shape, 
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The properties of the filters are clearly seen in the Fourier space. By embed¬ 
ding the eigenvector in the M x N matrix 
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we dehne the discrete Fourier transform as 
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and the completeness relations (perfect reproduction in the language of dig¬ 
ital filtering) 

^EIV$l 2 = l- (16) 

k =1 

The above properties Eqs. (fT5j) and (fT5j) are satisfied for arbitrary orthonor¬ 
mal vectors The SSA algorithm employs the eigenvectors of the eigen¬ 

value equation for X T X 

X T Xv {k) = X k v {k) , (17) 

and this criterion means that the quantities 

u (fc)T X T Xv {k) = \ k , (18) 

are maximized preserving the orthonormality of the vector [9]. In this 
sense, the SSA algorithm is the principal component analysis in the space 
of the lagged vectors. In the next subsection, we discuss the symmetry 
properties of the eigenvectors {rn fe '} and it’s implication to the 2D image 
decomposition. 

2.2 Symmetry properties of the filters 

We assign the vector components for the rectangular-shaped moving window 
with K — m x n in a natural order as in Eq. (1511 . We assume the periodic 
boundary condition for the image data. Then the lag-covariance matrix 
C = X r X/[L — K + 1) becomes to be the bisymmetric matrix, satisfying 
the relation 

Cij = CK-i+l,I<-j+l • (19) 

This can be seen as follows. The moving window in Eq. (JH]) slides from left to 
right and up to bottom sweeping the whole image data A. The component 
Cij of the matrix C is calculated by picking and summing up the matrix 
elements of A, corresponding to the position v t and Vj in the window, and 
this is the same for CK-i+i,K-j+i with Vx-i+i and Vx-j+i- The location of 
the components Vi and VK-i+i is just point symmetric about the center of 
the window. This also is the same for the relation between v 3 and vx-j+i- 
Thus the relative location between v, t and v 3 is the same as that between 
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Vk-i+i and vx-j+i- Because of this, the matrix element C t] is equal to 
CK-i+i,K-j+ 1 - This property holds for higher dimensional cases if we adopt 
the periodic boundary condition. For one dimensional case, a single-channel 
time series for instance, the matrix C is a symmetric Toeplitz matrix. 

Thus the matrix C commutes with the exchange matrix J 

[C, J]= 0 , (20) 

and if there is no degeneracy for the eigenvalues of C, the eigenvectors 
are either symmetric or antisymmetric 

Jv W = ±y (k) _ 


This means that the filters generated with the eigenvectors have cor¬ 

responding symmetry. In particular, for the case of square form m = n the 
matrices expressing the filter F ^ 
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are the centrosymmetric or skew-centrosymmetric matrices corresponding to 
the symmetric — v ( k ) or antisymmetric Jv^ k > = —v^ eigenvectors, 

respectively. 

As an example, we have carried out the SSA decomposition of the im¬ 
age data ’Building’ with rather small filter of square form 3x3. The nine 
eigenvalues of the matrix X T X are ordered from larger to smaller ones and 
are shown in Fig. 1. The five eigenvectors corresponding to the larger eigen¬ 
values are shown in Table 1. Here, the components of the eigenvectors are 
ordered in the 3x3 filter as in Eq. fl21|) . As seen, the eigenvectors r/ 1 ), 
rb 4 ) and v ^ are symmetric, while v ® , v ® are antisymmetric. 

To see the details of the decomposition characteristics, we have analyzed 
the filters as the differential operation. Let us assume that the image pixels 
are located on the lattice points ( ih,jk ), i , j G Z with mesh length h and 
k for both x and y directions, respectively. Here, the oxaxis is taken for the 
horizontal direction from left to right and the r/-axis for the vertical direction 
from up to bottom. The pixel value can be expressed with the function 
f{x = ih,y = jk ) on the lattice points, and the operation of the 3x3 filter 
F( k > can be expressed as 
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Figure 1: The nine eigenvalues for the SSA decomposition of the image 
’Building’ with 3x3 filters. The black circles(white triangles) are the eigen¬ 
values corresponding to the symmetric (antisymmetric) eigenvectors. 


v{ k) f{x -h,y + k)+ k) f(x - h, y ) 
vi k) f( x ~h,y~k) + v[ k) f(x, y + k ) 
vi k) f(x,y) + v { 6 k) f(x,y- k) 
v {k) f(x + h,y + k) + v$ k) f(x + h, y) 
v ( 9 k) f{x + h,y- k) . 

By carrying out the Taylor expansion we obtain 

f(x,y ) —> f(x,y)(J2 v i) 

i=l 

+ hf x (—V 1 -V 2 -V 3 +V 7 +V 8 +V 9 ) 

+ kf y (vi-v 3 +v 4 -v 6 + v 7 -v 9 ) 

+ h fxx CFl +^ 2+^3 + V 7 + Vg + Vg)/2 
+ hkf xy (-v ! + v 3 + v 7 -v 9 ) 

+ k 2 fyy(v 1 + l >3 + V 4 + V§ + ^7 + V g)/2 


f(x,y) — 
+ 
+ 
+ 
+ 


( 22 ) 
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Table 1: The larger five eigenvectors of the matrix X T X for the SSA de¬ 
composition of the photo image ’Building’ with 3x3 moving window. 


ih 1 ) 


v (3) 

F 4 > 


0.3327 

0.3747 

-0.4151 

-0.2671 

0.0899 

0.3338 

-0.0197 

-0.4336 

0.4223 

0.2913 

0.3326 

-0.4168 

-0.3732 

-0.2390 

0.3057 

0.3336 

0.4307 

-0.0198 

-0.1967 

-0.4716 

0.3348 

0.0000 

- 0.0000 

0.5558 

-0.4272 

0.3336 

-0.4307 

0.0198 

-0.1967 

-0.4716 

0.3326 

0.4168 

0.3732 

-0.2390 

0.3057 

0.3338 

0.0197 

0.4336 

0.4223 

0.2913 

0.3327 

-0.3747 

0.4151 

-0.2671 

0.0899 


+ higher order terms , (23) 

where the vector components are abbreviated as Vj and the f x = df(x, y)/dx 
etc. The coefficients for each terms are shown in Table 2 for dominant five 
filters. The filters F^ generated from the symmetric (antisymmetric) eigen¬ 
vectors can be seen as the derivative filters of the even (odd) order which is 
simply a consequence of the centrosymmetric (skew-centrosynnnetric) nature 
of the filters . 

Table 2: The coefficients appeared in Eq. (l23p for larger five eigenvectors. 
These clearly show how the filters work for the original image. 


Filter 

/ 

fx 

fy 

f XX 

fxy 

fyy 

1 

3.000 

0.000 

0.000 

0.999 

- 0.000 

0.999 

2 

0.000 

0.123 

2.444 

0.000 

0.000 

0.000 

3 

0.000 

2.444 

-0.123 

0.000 

0.000 

0.000 

4 

-0.005 

0.000 

0.000 

-0.084 

0.056 

-0.703 

5 

0.004 

0.000 

0.000 

0.687 

0.432 

-0.076 


The dominant filter F-0 works as smoother of the pixel values. The sec¬ 
ond and the third filters F^ and F^ are generated from the antisymmetric 
eigenvectors and work as the differential filter in y and x directions, respec¬ 
tively. The filters and F^ are generated from the symmetric eigen¬ 
vectors and the dominant components are the second derivative for various 
directions. 









( 5 ) ( 6 ) ( 7 ) ( 8 ) ( 9 ) 


Figure 2: The original image (upper left) and its SSA decomposition through 
3x3 filters for the photo image ’Building’. The images (k) are obtained by 
the two-step filtering with F( k \ 


The image decomposition for the ’Building’ are shown in Fig. 2. We can 
see that the results of the decomposition reflect the characteristics of each 
filter. For the decompositions (2) and (4), the horizontal lines are enhanced 
corresponding to the large f y ,f yy components as in Table 2. On the other 
hand, the vertical components are enhanced in the image (3) and (5). As 
seen, each filter has it characteristics and the SSA decompositions reflect 
their characteristic properties. The images (7) ~ (9) correspond to smaller 
eigenvalues and they exhibit somewhat featureless pattern due to the large 
high-frequency components. These nine images are added to reproduce the 
original one in the SSA algorithm. 

2.3 SSA decomposition of images with noise 

To see how the SSA algorithm works for the decomposition of the image data 
with noise, we adopt an example ’Lenna’ with 256 x 256 pixels. Each pixel 
takes the values 0 ~ 255. The Gaussian noise 13.3A r (0,1) is added for every 
pixels. We have carried out SSA decomposition with a square filter of the 
size 11 x 11. The larger fifteen eigenvalues of the matrix X T X are shown 
in Fig. 3, where the black circles (white triangles) represent the eigenvalues 
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corresponding to the symmetric (antisymmetric) eigenvectors. 



Figure 3: The larger fifteen eigenvalues of X T X for the SSA decomposition 
of the image ’Lenna’ with noise. The black circles (white triangles) are the 
eigenvalues with symmetric (antisymmetric) eigenvectors. 

Similar to the case of the ’Building’ image, we have examined how each 
filter works for the image by examining the coefficients in Eq. (T23T) extended 
for the case of 11 x 11 filter. In this case, since the filter size is large, we have 
replaced the mesh lengths h and k by h /3 and k/ 3, respectively as, 

f(x,y) —> v[ k) f{x-5h/3,y+5k/3)+v ( 2 k) f(x-5h/3,y+4:k/3)+ .(24) 

The coefficients /, hf x , kf y etc. are shown in Table 3. 

For the ’Lenna’ image with noise, the filter characteristics are quite similar 
to the case of the ’Building’ photo image. The first filter F ( - 1) works as the 
smoother for the original image. The smoothing operation works stronger 
than the case of ’Building’, since the filter size is larger 11 x 11 in this case. 
The second and the third Liters work as the differentiation to the x and y 
directions, respectively. For the fourth and the fifth Liters F ^ and F^\ the 
expansion coefficients f xx and f xy are large as seen in Table 3. 

For the Liters corresponding to the smaller eigenvalues, the derivative 
components are fragmented and the Lltered images gradually become to be 
structureless including the the higher frequency components. 
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Table 3: The coefficients appeared in Eq. (123 D for larger five eigenvectors for 
the ’Lenna’ with noise. In this case, since the filter size is large, we have 
replaced the mesh lengths h and k by h /3 and k /3 as in Eq. (12411 . The 
filters generated from the symmetric (antisymmetric) eigenvectors work as 
the derivative filters of the even (odd) order. 


Filter 

/ 

fx 

fv 

fxx 

fxy 

fyy 

1 

11.000 

0.000 

0.000 

6.079 

-0.021 

6.101 

2 

0.000 

11.337 

1.672 

0.000 

0.000 

0.000 

3 

0.000 

-1.669 

11.133 

0.000 

0.000 

0.000 

4 

0.067 

0.000 

0.000 

5.092 

2.654 

0.555 

5 

0.005 

0.000 

0.000 

-1.245 

11.140 

0.667 


In Fig. 4, we have shown the ’Lenna’ image with and without noise 
together with the dominant five SSA decomposition images with noise. Cor¬ 
responding to the filter characteristics as seen in the Table 3, we can see that 
the edge lines in various directions are enhanced and are extracted from the 
original image. For the images (2) and (3), the vertical and the horizontal 
lines are enhanced. The filters F*A and have large f xx and f xy compo¬ 
nents, and hence the finer details of the edges in the vertical or the diagonal 
lines are enhanced for the images (4) and (5). 

In order to see the relation between the SSA decomposition and the noise 
components, we show all the 121 eigenvalues for the ’Lenna’ image with and 
without the Gaussian noise in Figj5] The added noise enhances the smaller 
eigenvalues coming from the high-frequency components. 

We have also calculated the RMS distance D(t) between the original 
noiseless image a tJ and the partial sum of the SSA decomposition [?], 


m 


\ 


i,j k =1 


(25) 


The results are shown in Fig. [6] Obviously, for the noiseless case, the distance 
converges to zero as £ —» 121, which merely means that the sum of the SSA 
decomposition reproduce the original image. For the case with noise term, 
the RMS distance decreases as i : 0 —* 37, while it increases as £ : 37 —> 121. 
The RMS distance takes minimum value at l = 37. We could discard the 
components {d^] k >37} in order to obtain better images. The image 
corresponding to the partial sum k = 1 ~ 37 is also shown in Fig. [TJ The 
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Lenna.bmp with Noise filterl ~37 



Figure 4: Upper left and central images are the ’Lenna’ with and without 
Gaussian noise. The image is decomposed into 11 x 11 components with SSA 
algorithm. The upper right image is the partial sum of the SSA decomposi¬ 
tion k — 1 ~ 37. The lower images are the dominant five decompositions. 

RMS distance between noiseless and noisy images is D{1 = 121) = 13.3 which 
is just the standard deviation of the added noise. The RMS distances for the 
partial sum of the SSA decomposition k = 1 ~ 37 with and without noise are 
D(( 37)| w jth noise 8.42 and D(( 37)| no iseiess 5.21, respectively. 

So far, much efforts have been devoted to the image restoration. Various 
linear or non-linear filtering operations have been examined to improve the 
noisy images [1,5]. The performance of the filtering operation depends on the 
properties of the image itself or the nature of the noise. In the SSA, the noisy 
images are exactly decomposed into arbitrary number of components and 
each image decomposition has its own characteristics. Thus, by combining 
the conventional nonlinear denoising algorithms and SSA decomposition, we 
expect to develop better denoising algorithms. Since the main purpose of 
this paper is to show the characteristic feature of the SSA decomposition of 
the 2D image data, we do not pursue these points further. 
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Figure 5: The 11 x 11 eigenvalues are shown for the SSA decomposition of 
the image ’Lenna’. The dashed and the solid lines are the results with and 
without the Gaussian noise. 

3 Summary and Conclusions 

The filtering interpretation of the SSA algorithm enabled us to study the 
characteristics of the image data decomposition in detail. The periodic 
boundary condition and the rectangular moving window leads to the bisym- 
metric lag-covariance matrix C and its eigenvectors are either symmetric or 
antisymmetric. For the symmetric (antisymmetric) eigenvectors, the corre¬ 
sponding filter F^ are centrosymmetric (skew-centrosynnnetric) matrices 
and they are interpreted as the derivative filter of even (odd) order. The 
SSA algorithm can be understood as the optimal and adaptive generation of 
the filters F ® and the two-step point-symmetric operation of these filters 
to the reference points of the original image. The two-step filtering ensures 
the reality of the filters in the Fourier space and it introduces no distortion 
or shift to the data. Also the completeness of the eigenvectors ensures the 
perfect reproductivity of the SSA decomposition. Each image components 
in the SSA decomposition exhibit characteristic features, which comes from 
the properties of the adaptive filters generated from the eigenvectors of the 
lag-covariance matrix. We briefly discussed the relation of the SSA decompo¬ 
sition and the noise reduction of the images. Much efforts have been devoted 
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Figure 6: The RMS distance D(£) between the original image and the partial 
sum of the SSA decomposition k — 1 ~ i as in Eq. (125|) . The solid and the 
dashed lines are the results with and without Gaussian noise 13.6iV(0,1). 
The horizontal line represents the RMS distance between the noiseless and 
noisy images. 


so far for the image denoising. We expect that the filtering approach for SSA 
algorithm enables us to develop an novel denoising algorithm by combining 
the SSA decomposition with conventional nonlinear filtering algorithms. The 
study in this direction is now under way. 
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